function test
 a=[1/2 1/3 1/4];
b=[1/5 1/6 1/3];
c=cross(a,b);
disp("a x b=");
disp(c);

A=pauliVectorProduct(a)/2;
B=pauliVectorProduct(b)/2;
disp("A=");
disp(A);
disp("B=");
disp(B);
disp("iA,iB in su(2)?");

C1=expm(-commutatorMatrix(A,B));
disp("[A,B]=");
disp(commutatorMatrix(A,B));
C2=u_rot(c);

disp("e^-[A,B]=");
disp(C1);

disp("u(a x b)=");
disp(C2);
disp("sqrtm of u(a x b)=");
disp(my_sqrt(C2));
disp("or");
disp(sqrtm(C2));

endfunction

function result=commutatorMatrix(A,B)
 if(size(A)~=size(B))
  error('Incompatible matrix size');
 endif 
 result=A*B-B*A; 
endfunction

